// polyselector.cc
// Program for selecting polynomials for the General Number
// Field Sieve (GNFS).
// It creates an irreducible  poynomial of wanted degree....
// It is based on the NTL library, which is based on GnuMP

#include <iostream>
#include <ctime>
#include <string>
#include <NTL/ZZ.h>
#include <NTL/RR.h>
#include <NTL/ZZX.h>
#include <NTL/ZZ_pXFactoring.h>
#include "nfsinput.h"

NTL::RR DickmanRho(NTL::RR u) {
	// Calculation of the Dickman rho functions,
	// which is the probability that an integer
	// has some greatest prime factor as F(x) reaches
	// some given bound!
	if(u < 1.0) {
		return NTL::to_RR(1.0);
	}

	NTL::RR K = NTL::ceil(u);
	NTL::RR k,tmpZeta;
	NTL::RR zeta = k-u;
	NTL::RR tmpc;
	int tableSize=55;
	NTL::RR* c = new NTL::RR[tableSize];
	NTL::RR* Newc = new NTL::RR[tableSize];

	if(K>=NTL::to_RR(2)) {
		c[0] =  1.0 - NTL::log(NTL::to_RR(2.0));

		for(int i=1; i<tableSize;i++) {
			tmpc = NTL::to_RR(i)*NTL::power(NTL::to_RR(2.0),i);
			c[i] = NTL::inv(tmpc);
		}
		k=NTL::to_RR(2);
	}
	else {
		c[0] =  1.0;

		for(int i=1; i<tableSize;i++) {
			c[i] = NTL::to_RR(0);
		}
		k=NTL::to_RR(1);
	}
	
	while (k < K) {
		for (int i=1; i<tableSize; i++) {
			tmpc = NTL::to_RR(0.0);
			for (int j=0; j<i; j++) {
				tmpc += c[j]/(NTL::to_RR(i)*NTL::power(k+1,i-j));
			}
			Newc[i] = tmpc;
		}
		tmpc = NTL::to_RR(0.0);
		for (int j=1; j<tableSize; j++)
			tmpc += Newc[j]/NTL::to_RR(j+1);
		Newc[0] = tmpc/k;
		for(int g=0;g<tableSize;g++) {
			c[g]= Newc[g];
		}
		k++;
	}
	tmpc = NTL::to_RR(0.0);
	tmpZeta=NTL::to_RR(1.0);
	for (int i=0; i<tableSize; i++) {
		tmpc += c[i]*tmpZeta;
		tmpZeta *= zeta;
	}
	delete[] c;
	delete[] Newc;
	return tmpc;
}

NTL::RR CalcAlpha(NTL::ZZX f,long sampleBound, NTL::ZZ* primeBase, int primeBound ) {
	// ** It calculates the alpha function denoted in  
	// ** Murphy's thesis.
	
	NTL::ZZ p,j,a,b; 
	NTL::ZZ disc_f=NTL::discriminant(f);
	NTL::ZZ a_d = NTL::LeadCoeff(f);
	NTL::ZZ q_p;
	NTL::RR alpha=NTL::to_RR(0);
	NTL::RR cont_p;
	NTL::vec_ZZ_p roots;
	NTL::ZZ realBound= NTL::to_ZZ(NTL::pow(NTL::to_RR(sampleBound),NTL::to_RR(2)));

	for(int i=0;i<primeBound;i++) {
		p=primeBase[i];
		q_p=NTL::to_ZZ(0);	
		if(disc_f%p != 0) {

			for(int i=0;i<sampleBound;i++) {
				for(;;){
					a=NTL::RandomBnd(realBound);
					if (a!=NTL::to_ZZ(0))
						break;
				}	
				b=NTL::RandomBnd(realBound);
				while(NTL::GCD(a,b)!= 1) {
					b=NTL::RandomBnd(realBound);
				}
				
			 	j= EvalF(f,a,b);
			 	while (j%p==0) {
					j/=p;
					q_p++;
				}
			}
			cont_p=NTL::MakeRR(q_p,0)/NTL::MakeRR(NTL::to_ZZ(sampleBound),0);
		}
		else {
			// Contribution from 
			// the real roots...
			for ( j=NTL::to_ZZ(0); j<p ; j++) {
				if( EvalfmodP(f,j,p)==0 ){
					q_p++;
				}
			}
			// Contribution from 
			// the projective roots...
			if(a_d%p == 0) {
				q_p++;
			}
			cont_p = NTL::MakeRR(q_p,0)/(NTL::MakeRR(p,0)-1);
		}
		NTL::RR tmpP = NTL::MakeRR(p,0);		
		alpha += (NTL::inv(tmpP-1) - cont_p)*NTL::log(tmpP);
	}							      
	return alpha;
}

NTL::RR RatePolynomial(NTL::ZZX f,long sampleBound, NTL::ZZ* primeBase, int primeBound ){
	// compute the rating of a polynomial according to
	// equation 5.6 in Murphy's thesis
	int K = 1000;
	NTL::RR logBound = NTL::log(NTL::to_RR(sampleBound));
	NTL::RR rating = NTL::to_RR(0.0);
	NTL::RR pi = NTL::ComputePi_RR();
	NTL::RR u;
	NTL::RR alpha = CalcAlpha(f,sampleBound,primeBase,primeBound);
	NTL::RR theta = pi/NTL::to_RR(2*K);
	
	for (int i=0; i<K; i++) {
        	u = (NTL::log(NTL::fabs(EvalF(f,NTL::cos(theta),NTL::sin(theta)))) + alpha)/logBound;
		rating += DickmanRho(u);
		theta += pi/NTL::to_RR(K);
	}
	return rating;
}

NTL::ZZX MakeBase_mPolynomial(NTL::ZZ n,NTL::ZZ m,int degree){
	// ** It returns a polynomial, created by a base-m 
	// ** expansion of the number.
	// ** n: the number to be expanded
	// ** m: is the basenumber
	// ** degree: the degree of the resulting polynomial
	NTL::ZZ newN=n;
	NTL::ZZ tmp,coeff;
	NTL::ZZX res;
	
	// set size of polynomial
	res.rep.SetLength(degree+1);
	
	// make the Base-m representation
	for(int i=degree; i>0 ; i--) {
		tmp = NTL::power(m,i);	
		coeff = newN/tmp;
		res.rep[i] = coeff;
		newN = newN - (coeff*tmp);
	}
	res.rep[0] = newN;
	NTL::ZZ m_2 = m/NTL::to_ZZ(2);
	
	for(int i=degree; i> 1; i--) {
		if(res.rep[i] > m_2) {
			res.rep[i] -= m;
        		res.rep[i-1] += 1;
		}
	}
	return res;
}

void GenPrimes(NTL::ZZ* primes, int n) {
	// Generate the first n primes and
	// store them in 'primes'
	NTL::PrimeSeq p;
      	primes[0]= NTL::to_ZZ(2);
	
	for(int i=1;i<n;i++) {
      		primes[i]= p.next();
	}
}


NTL::ZZX MakePolynomial(NTL::ZZ n, NTL::ZZ cofactor, int primebound,char* polyFile) {
	// ** Creates one suitable polynomial for the GNFS

	// Find a suitable degree of the polynomial
	int d = (int) floor(pow(((double) NTL::NumBits(n))/3.322,1.0/3.0)); 
	if(d<3) {
		d=3;
	}
	NTL::ZZ m,a_dmax,a_dmin,a_d,a;
	NTL::ZZX test;
	NTL::RR d_inv = 1/(NTL::MakeRR(NTL::to_ZZ(d),0));
	nfsPoly best=readPoly(polyFile);
	nfsPoly murphy=readPoly("murphy.poly");
	nfsPoly candidate;
	
	// construct array of primes below primeBound
	// Using approx of number of primes ~x/ln(x)
	int numPrimes = (int) floor((double) primebound/log((double) primebound));
	NTL::ZZ* primes = new NTL::ZZ[numPrimes];
	GenPrimes(primes, numPrimes);
	
	// estimate a value for m ( m=n^(1/d)/2 )  
	NTL::RR l = 1/(NTL::MakeRR(NTL::to_ZZ(d),0)+1);
	m=NTL::FloorToZZ( NTL::pow( NTL::MakeRR(n,0),l));
	
	// pick an interval for a_d
	a_dmin =  m/10000;
	a_dmax =  m/2;
	
	// find combinations of a_d and m which gives
	// small a_{d-1}
	a_d = a_dmin + (a_dmax-a_dmin)/2;	
	
	NTL::SetSeed(NTL::to_ZZ(time(NULL)));
		
	
	while(1) {
		a = NTL::RandomBnd(a_d);
		if (a==0) {
			a=1;
		}
		// use the given cofactor and make sure that a_d has
		// the given cofactor!
		if (cofactor!=0 && (a- a%cofactor)!=0) {
			a = a - (a%cofactor);
		}
		candidate.m = NTL::FloorToZZ(NTL::pow(NTL::MakeRR(n,0)/NTL::MakeRR(a,0),d_inv));
		candidate.f = MakeBase_mPolynomial(n,candidate.m,d);
		candidate.alpha = CalcAlpha(candidate.f,1000,primes,numPrimes);
		// if this polynomial is better than the best one yet
		// keep it!!!
		if (best.alpha == 0 || candidate.alpha < best.alpha) {
			std::cout << "* * * * * * * * * * * * *" << std::endl;
			std::cout << "* * *  New BEST!!!  * * *" << std::endl;
			std::cout << "* alpha: " << candidate.alpha << std::endl;
			std::cout << "* * * * * * * * * * * * *" << std::endl;
			best.f=candidate.f;
			best.m=candidate.m;
			best.alpha=candidate.alpha;
			best.rating=candidate.rating;
			storePoly(best,polyFile);
			//printPoly(best);
		}
	}	

	// clean up...
	delete[] primes;
	return best.f;
}


int main (int count, char** args){
	NFS_instance prob = readData(args[1]);

	std::string polyfile=args[1];
	polyfile+=".poly";
	NTL::RR::SetPrecision(NTL::NumBits(prob.n));
		
	std::cout << MakePolynomial(prob.n,prob.cofactor,NTL::to_long(prob.primebound),(char*) polyfile.c_str()) << std::endl;
}

